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We investigate the solvent-accessible area method by means of Metropolis simulations of the brain 
peptide Met-Enkephalin at 300 if. For the energy function ECEPP/2 nine atomic solvation param- 
eter (ASP) sets are studied. The simulations are compared with one another, with simulations with 
a distance dependent electrostatic permittivity e(r), and with vacuum simulations (e = 2). Par- 
allel tempering and the biased Metropolis techniques RMi are employed and their performance is 
evaluated. The measured observables include energy and dihedral probability densities (pds), inte- 
grated autocorrelation times, and acceptance rates. Two of the ASP sets turn out to be unsuitable 
for these simulations. For all other systems selected configurations are minimized in search of the 
global energy minima, which are found for the vacuum and the e(r) system, but for none of the 
ASP models. Other observables show a remarkable dependence on the ASPs. In particular, we find 
three ASP sets for which the autocorrelations at 300 K are considerably smaller than for vacuum 
simulations. 

PACS numbers: 05.10.Ln, 87.16-v, 87.14.Ee. 



I. INTRODUCTION 

In nature biomolecules exist in the environment of 
solvents, thus the molecule-solvent interactions must be 
taken into account. It is very computer time consuming 
to simulate models for which the molecules of the sur- 
rounding water are treated explicitly. Therefore, a num- 
ber of approximations of solvent effects have been devel- 
oped. In the solvent-accessible area approach 0,0,0] it is 
assumed that the protein-solvent interaction is given by 
the sum of the surface area of each atomic group times the 
atomic solvation parameter (ASP). The choice of a set of 
ASPs (also called hydrophobicity parameters or simply 
hydrophobicities) defines a model of solvation. However, 
there is no agreement on how to determine the univer- 
sally best set of ASPs, or a at least the best set for some 
limited purpose. For instance, eight sets were reviewed 
and studied by Juffer et al. [4( and it was found that 
they give rather distinct contributions to the free energy 
of proteins folding. 

In this paper we investigate how different ASP sets 
modify the Metropolis simulations of the small brain pep- 
tide Met-Enkephalin (Tyr-Gly-Gly-Phe-Met) at 300 K. 
The reason for the choice of Met-Enkephalin is that its 
vacuum properties define a reference system for testing 
numerical methods, e.g. [E El EJ- Therefore, 
Met-Enkephalin appears to be well suited to set refer- 
ences for the inclusion of solvent effects as well, but we 
are only aware of few articles [TTL Il2| , which comment on 
the modifications due to including a solvent model. 

We set our simulation temperature to 300 K, because 
room temperature is the physical temperature at which 
biological activity takes place. Most of the previous sim- 



ulations of Met-Enkephalin in vacuum were performed 
at much lower temperatures or employed elaborate mini- 
mization techniques with the aim to determine the global 
energy minimum (GEM) . Only recently [T(il ] it was shown 
that the GEM is well accessible by local minimization 
of properly selected configurations from an equilibrium 
time series at 300 K. Precisely this should be the case for 
a GEM which is of relevance at physical temperatures. 



For our simulations we use the program package 
SMMP El (Simple Molecular Mechanics for Proteins) 
together with parallel t emp ering 

El 13 M ( PT ) and 
the recently introduced 10] biased Metropolis technique 
RMi (rugged Metropolis - approximation 1). SMMP 
implements a number of all-atom energy functions, de- 
scribin g th e intramolecular interactions, and nine ASP 
sets Ulnl El IH IH Him HI to model the molecule 
solvent interactions. We use the ECEPP/2 [24[ (Empir- 
ical Conformational Energy Program for Peptides) en- 
ergy function with fully variable lu angles and simulate 
all nine ASP sets. For comparison we simulate also Met- 
Enkephalin in vacuum and with the distance dependent 
electrostatic permittivity e(r) of Ref. 25]. 



The paper is organized as follows: The energy func- 
tions and Metropolis methods used are explained in 
Sec. [n] In Sec. IIIII we present our results from simu- 
lations of the brain peptide Met-Enkephalin. Summary 
and conclusions are given in Sec. II VI 
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II. MODELS AND METHODS 
A. ASP sets 

In all-atom models of biomolecules the total confor- 
mational energy of the intramolecular interactions E\ is 
given as the sum of the electrostatic, the Lennard-Jones 
(Van der Waals), the hydrogen bond, and the torsional 
contributions, 
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Here r,j is the distance between atoms i and j, qi and 
qj are the partial charges on the atoms i and j, e is the 
electric permittivity of the environment, Aij, Bij, Cij 
and Dij are parameters that define the well depth and 
width for a given Lennard-Jones or hydrogen bond inter- 
action, and 4>k is the fcth torsion angle. The units are 
as follows: distances are in A, charges are in units of the 
electronic charge and energies are in kcal/mol. 

One of the simplest ways to include interactions with 
water is to assume a distance dependent electrostatic per- 
mittivity according to the formula [2^, [2|J 



e(r) = D 



D-2 



[(sr) 2 + 2sr + 2] 



(2) 



Empirical values for the parameters D and s are cho- 
sen, so that for large distances the permittivity takes the 
value of bulk water, e = 80, and the value e = 2 for 
short distances, i.e. for the interior of the molecule. Ap- 
proximating solvation effects in this way is implemented 
as an option in SMMP. It allows to include solvation ef- 
fects without any significant slowing down over the vac- 
uum simulation with e = 2. The approach is clearly an 
oversimplification, because atoms which are close to each 
other do not necessarily have to be simultaneously in the 
interior of the molecule. Reversely, two atoms which are 
separated by a large distance may still be in the interior 
of the molecule. More elaborated approaches are asked 
for. 

If the molecule-solvent interaction is proportional to 
the surface area of the atomic groups, it is given by the 
sum of contributions of a product of the surface area 
of each atomic group and the atomic solvation parame- 
ter H, 



E, 



sol 



(3) 



Here E so \ is the solvation energy and the sum is over all 
atomic groups. A$ is the solvent accessible surface area 
and <7j the atomic solvation parameter of group i. The 
choice of a set of ASPs a* defines a model of solvation. 
There are nine sets of ASPs in the SMMP package, which 



TABLE I: Atomic solvation parameter sets implemented in 
SMMP. The first column gives the value of the SMMP pa- 
rameter itysol and the second column the letter code used in 
SMMP. In the author column we give also the year of publica- 
tion. The last column indicates the method used as explained 
in the text. 







Authors 




1 


OONS 


Ooi, Obatake, Nemethy, Scheraga [17] 1987 


v/w 


2 


JRF 


Vila, Williams, Vasquez, Scheraga [20] 1991 


v/w 


3 


WE92 


Wesson and Eisenberg [21] 1992 


v/w 


4 


EM86 


Eisenberg and McLachlan [3J 1986 


o/w 


5 


SCH1 


Eisenberg, Wesson, Yamashita [181 1989 


o/w 


6 


SCH2 


Kim [19], see [4] 1990 


o/w 


7 


SCH3 


Wesson and Eisenberg [21] 1992 


v/w 


8 


SCH4 


Schiffer, Caldwell, Kollman, Stroud [22] 1993 


v/ws 


9 


BM 


Freyberg, Richmond, Brown [23] 1993 


cla 



we list in Table||J Columns one and two of this table gives 
the notations used in SMMP to identify the different sets. 

Eisenberg and McLachlan were the first to deter- 
mine a set of ASPs (itysol= 4, EM86 in the SMMP 
notation). For this, they considered the process of trans- 
ferring atoms or groups of atoms from the interior of a 
protein to aqueous solution and used transfer energies of 
amino acids from n-octanol to water as reported in [27j . 
The ASPs are then determined by least-square fitting. 
Octanol is chosen, because it apparently resembles the 
interior of a protein. With the exception of [2^] and [2^] , 
the other authors used similar methods with the major 
variation that instead of transfer energies with respect to 
octanol- water (o/w) also transfer energies with respect 
to vacuum- water (v/w) were used (for early determina- 
tion of v/w transfer energies see [28|, 29]). The last col- 
umn of Table [I] indicates whether the transfer energy is 
o/w or v/w. In the chronological order [3L ITsL l2ll | Eisen- 
berg and collaborators contributed the parameter sets 
EM86, SCHLWE92 and SCH3. Scheraga and collab- 
orators 0, |2(J contributed the parameter sets OONS 
and JRF. Here it should be noted that some of the 
original ASP sets were modified in course of time. For 
itysol= 1, ... ,8 SMMP implements the parameters as 
reviewed and tabulated in Ref. 3 , where in turn the sets 
SCH1 to SCH4 are simply taken from Schiffer et al.JH. 
Table 1 of SMMP lists the implemented ASPs for 
itysol = 1, . . . , 8. 

Somewhat special cases are the ASP sets SCH4 22] 
and BM [23| • SCH4 was determined by comparison of the 
crystal structure from molecular dynamics simulations of 
small peptides and proteins in explicit water with similar 
simulations using an ASP solvation term (v/ws). The 
BM set of SMMP relies on a specific classification (cla) 
of atomic groups, where for all non-hydrogen atoms the 
solvation coefficients are set to 1 kcal/mol per A 2 . 



3 



B. Metropolis methods 

For the updating of our systems we use PT with two 
processors, one running at 300 K and the other at 400 K. 
This builds on the experience |10| | with vacuum simula- 
tions of Met-Enkephalin for which the following observa- 
tions are made: 

1. The integrated autocorrelation time r; n t (defined 
below in this section) increases from 400 K to 300 K 
by a factor of ten for the (internal) energy and by 
factors of more than twenty for certain dihedral an- 
gles. 

2. The energy probability densities (pds) at 300 K and 
400 K overlap sufficiently, so that the PT method 
works, leading to an improvement factor of about 
2.5 in the real time needed for the simulation (see 
Table I of Ref. [13). 

A brief description of the PT algorithm is given in the 
following. PT performs n canonical MC simulations at 
different /3-values with Boltzmann weight factors 
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where Po < Pi < ■ ■ ■ < Pn-2 < ftn-i arid a configuration 
is denoted by k. PT allows the exchange of neighboring 
/3-values 

Pi-t <-» Pi for i = 1, . . . , n - 1 . (5) 
These transitions lead to the change 

-AH = -Pi^ (E^ - Efl) - P t (e^\ - Ef) 
= {Pi -P^) (e^-E^ (6) 



which is accepted or rejected according to the Metropolis 
algorithm, i.e. with probability one for AH < and with 
probability exp(-Aff) for AH > 0. 

For the vacuum system the performance of the PT 
simulation was improved by an additional factor of two 
in Ref.jlOJ by using a first approximation, called RMi, 
to the rugged Metropolis scheme introduced there. A 
(short) simulation at 400 K was used to obtain estimates 
~Pj(vj), j = 1, . . . , 24 of the pds of the 24 dihedral angles, 
which were then fed into the simulation. For a configu- 
ration change k — > k' at temperature Tj the new config- 
uration accepted with the probability 
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in the RMi updating scheme. In the present paper we 
report the improvement due to this biased updating for 
some of the ASP sets. 

In the vacuum simulation it is possible to determine 
the GEM by minimizing selected configurations of the 




100 200 300 400 500 600 
t [32 sweeps] 

FIG. 1: Selection of configurations for local minimization from 
the energy time series at 300 K. The lower and upper straight 
lines indicate the quantiles Eq.i and Eo.g, respectively. 



300 K time series. Here we apply the same procedure 
to our PT simulation of the ASP sets introduced in the 
previous subsection: 

1. We determine the lower 10% quantile and the 
upper 10% quantile E , 9 of the energy distribution 
of our time series. This is done by sorting all ener- 
gies in increasing order and finding the values which 
cut out the lower and upper 10% of the data. For 
the statistical concepts see, e.g., Ref. [30| . 

2. We partition the time series into bunches of config- 
urations. A bunch contains the configurations from 
one crossing of the upper quantile Eq q to the next 
so that at least on crossing of the lower quantile 
Eq.i is located between the two crossings of So. 9. 
For each bunch we pick then its configuration of 
lowest energy. The idea behind this procedure is to 
pick minima of the time series, which are to a large 
degree statistically independent. In Fig. ^ the ar- 
rows indicate the energy values picked in that way 
from the first 600 configurations recorded in the 
RMi simulation of Ref. 10]. 

3. We run a conjugate gradient minimizer on all the 
selected configurations and thus obtain a set of con- 
figurations which are local energy minima. For the 
vacuum simulation [10( about 5% to 6% of the thus 
minimized configurations agree with the GEM. 

To determine the speed at which the systems equili- 
brate, we measure the integrated autocorrelation time 
Tint for the energy and each dihedral angle. In particular 
the integrated autocorrelation times are directly propor- 
tional to the computer run times needed to achieve the 
same statistical accuracy for each system. They thus de- 
termine the relative performance of distinct algorithms. 
For an observable / the autocorrelations are 



C(t) = (f Q f t )-(f)> 



(8) 
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where t labels the computer time. Defining c(t) — 
C(t)/C(0), the time-dependent integrated autocorrela- 
tion time is given by 



71* (t) = l + 2^e(t') 



(9) 



Formally the integrated autocorrelation time Ti n t is de- 
fined by ri n t = lim^oo T lnt (t). Numerically, however, this 
limit cannot be reached as the noise of the estimator in- 
creases faster than the signal. Nevertheless, one can cal- 
culate reliable estimates by reaching a window of t values 
for which T; n t(i) becomes flat, while its error bars are still 
reasonably small. This is the method we employ in the 
next section, see Ref. jU for a more detailed discussion 
of the integrated autocorrelation time. 



III. RESULTS 
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FIG. 2: The time-dependent integrated autocorrelation time 
for the energy at 300 K from our simulations of the vacuum 
system and two of the solvent models of Table Q 



A. Autocorrelations 

The PT simulations with temperatures To = 400 K 
and Ti = 300 K are performed on the system in vacuum 
(e = 2), with e(r) given by Eq. @ and for the nine 
ASP sets of Tabled The dihedral angles updated in our 
simulations are fully variable in the range from — n to n. 
We keep a time series of 2 16 = 65536 configurations for 
each replica (i.e., each of the two processor), in which 
subsequent configurations are separated by 32 sweeps. A 
sweep is defined by updating each dihedral angle once 
sequentially. Before starting with measurements 2 18 = 
262144 sweeps are performed for reaching equilibrium. 
Thus, the entire simulation at one temperature relies on 
2 2i + 2 is = 2,359,296 sweeps. On the Cray T3E, this 
takes about 14 hours for the vacuum system and 5 x 14 
hours for each ASP set. 

Results of the average energy, acceptance rates and in- 
tegrated autocorrelations times for the energy E variable 
are shown in Table ITD For the vacuum simulations and 
the ASP sets OONS and EM86 the time-dependent inte- 
grated autocorrelations times © are shown in Fig. [21 In 
each case a window of t values is reached for which n nt (t) 
does no longer increase within its statistical errors. In the 
case of the vacuum simulations it even does decrease, but 
this is not significant due to the statistical error. These 
windows are then used to estimate the asymptotic Tj nt 
values of Table ^ With the exception of the ASP sets 
JRF and BM, the integrated autocorrelations times of all 
other sets are determined in the same way. 

From Table |n] we see that the acceptance rates of the 
solvent models JRF and BM are much lower than for the 
other models. In essence the simulations of these two 
models get stuck, which implies that their integrated au- 
tocorrelation times cannot be measured. This is illus- 
trated in Fig. |3 for the time-dependent integrated au- 
tocorrelation time of the energy at 300 K. The function 
Tint (t) increases rapidly until it gets lost in the noise. The 



2000 




t [32 sweeps] 

FIG. 3: The time-dependent integrated autocorrelation time 
for the energy at 300 K from our simulations of the solvent 
models JRF and BM. 



pds of the dihedral angles of these two models are also 
erratic and the conclusion is that they cannot be used to 
describe Met- Enkephalin in solvent. 

The energy couples to all dihedral angles and its inte- 
grated autocorrelation time is characteristic for the en- 
tire system, while the integrated autocorrelation times 
of the single dihedral angles vary heavily from angle to 
angle. For all our systems, but JRF and BM, we show 
in Fig. 2] the integrated autocorrelation times at 300 K 
for the energy and all dihedral angles. The notation m, 
i = 0,1,..., 24 is used, where vq is stands in for the 
energy E and for i — 1, . . . , 24 the i?j are the dihedral 
angles as used in the SMMP computer program. The 
relationship of the Vi angles to the conventional nota- 
tion for dihedral angles and their residues is summarized 
in Table ITTT1 where it should be noted that the SMMP 
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TABLE II: Average energies (E) (Kcal/mol), acceptance rates and integrated autocorrelations times Tint for the energy are 
shown for simulations in vacuum (VAC), with e(r) of Eq. © and with the nine ASP introduced in Table [I] 





T = 400 K 


T — 300 K 




(E) 


acpt 


Tint 


(E) 


acpt 


Tint 


VAC 


7.07(03) 


0.167 


3.67(20) 


1.29(06) 


0.119 


19.9(1.6) 


e(r) 


-12.00(03) 


0.171 


2.92(10) 


-17.61(06) 


0.121 


14.35(75) 


OOJNb 


-13.80(01) 


0.195 


1.25(02) 


-17.70(02) 


0.143 


2.64(14) 


TTD 1 / 

Jrlr 


-311.69(44) 


0.058 




-319.08(40) 


0.046 




WE92 


-15.76(02) 


0.199 


1.30(03) 


-19.75(02) 


0.145 


2.94(07) 


EM86 


13.49(03) 


0.158 


4.71(21) 


8.03(06) 


0.116 


25.0(2.9) 


SCH1 


10.45(03) 


0.165 


3.72(22) 


4.95(06) 


0.119 


23.2(2.2) 


SCH2 


-18.33(02) 


0.212 


1.11(01) 


-21.83(01) 


0.160 


1.89(05) 


SCH3 


13.33(03) 


0.151 


4.59(34) 


8.33(06) 


0.112 


26.4(3.3) 


SCH4 


13.38(03) 


0.158 


4.35(17) 


7.85(05) 


0.115 


25.1(2.1) 


BM 


630.4(3.9) 


0.043 




610.6(3.0) 


0.037 






FIG. 4: Integrated autocorrelation times for the energies 
(«o = E) and the dihedral angles Vi, i = 1, . . . , 24 at T=300 K. 
The up-to-down order of the curves agrees at i = 10 with the 
order in the figure legend. 



notation differs from other literature 0, Q . 

In Fig.0|we see that for each dihedral angle Vi the inte- 
grated autocorrelation times Ti nt [wi] for the three solvent 
models OONS, WE92 and SCH2 are smaller than for the 
remaining systems, including the vacuum system. For 
the integrated autocorrelation time of the energy Ti nt [E] 
this observation is already obvious from Tabled In par- 
ticular, this means that the OONS, WE92 and SCH2 
models require far less statistics than the vacuum run 
for achieving the same accuracy of results. Using the 
TintfE 1 ] results of Table ITU we find a factor in the range 7 
to 10, which more than offsets the factor of 5 by which 
the ASP model simulations are slower than the vacuum 
simulation. In the following the solvation models OONS, 
WE92 and SCH2 define the "fast class" , while the other 
models shown in Fig. 21 constitute the "slow class" (the 
models JRF and BM are omitted from this classifica- 
tion). "Good" behavior of the models OONS and WE92 



has been previously observed |32| . 

The autocorrelation times in the fast class are so small 
that the resolution of 32 sweeps in our recorded time se- 
ries becomes too crude. Namely, autocorrelations over 
less than 32 sweeps are then not measured and the in- 
tegrated autocorrelation time approaches one as soon 
as autocorrelations stay within the range of 32 sweeps. 
To investigate this point further, we performed for the 
OONS, WE92 and SCH2 models simulations for which 
the configurations were recorded every four sweeps and 
the total statistics was reduced by the factor 1/8. In the 
new units of four sweeps the integrated autocorrelation 
time is larger by a factor which is bounded by 8 = 32/4. 
The bound is assumed, if there is no improvement due 
to integrating additional small fluctuations out (i.e. due 
to the additional configuration in between the 32 sweeps, 
which are now kept in the time series). 

For WE92 we report in the PT column of Table |TTI| 
the integrated autocorrelation times from the simulation 
with reduced statistics. For many dihedral angles the in- 
crease lies well below the factor of eight, showing that 
we gain in accuracy by averaging over small fluctuation 
within the range of 32 sweeps. On the other hand, noth- 
ing is gained by this extra averaging for several angles 
with large autocorrelations. In those case the simulations 
yield within their statistical errors the upper bound 8. 

To supplement the vacuum results of Ref. 01 > we re ~ 
peated the WE92 PT simulations by using estimates of 
the dihedral pds from 400 K as input for the biased up- 
dating of Eq. Q . These results are reported in the PT- 
RMi column of Table II I II As in the case of the vacuum 
simulations, we find an improvement of the PT perfor- 
mance by a factor of approximately two, which is also 
obtained for the other models of the fast class. For the 
slow class we checked by direct improvement of the origi- 
nal simulations on the ASP models EM86 and SCH4 and 
find again an acceleration by a factor of about two when 
we are using RMi updating. 
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TABLE III: Definitions of the dihedral angles together with 
their integrated autocorrelation times Ti nt at 300 K for simu- 
lations of WE92 with a statistics reduced by 1/8 and con- 
figurations recorded every four sweeps. PT denotes the 
400 K-300 K parallel tempering simulation. For PT-RMi the 
PT simulation is supplemented by the RMi biased updating 
Q with input pds from 400 K. The factor of the last column 
denotes the increase of the PT Ti nt over its values for the full 
WE92 simulations where configurations were recorded every 
32 sweeps (8 is the upper bound for this factor). 



var 


angle 


res [5, 8] 


res [131 


PT-RMi 


PT 


factor 


Vl 


x 1 


Tyr-1 


Tyr-1 


6.9 (1.1) 


11.6 (1.6) 


6.1 (0.9) 


«2 


x 2 


Tyr-1 


Tyr-1 


2.0 (0.2) 


3.1 (0.5) 


2.7 (0.5) 


V 3 


x 6 


Tyr-1 


Tyr-1 


1.0 (0.1) 


1.3 (0.2) 


1.3 (0.2) 


v 4 


<t> 


Tyr-1 


Tyr-1 


2.1 (0.2) 


2.6 (0.4) 


2.4 (0.4) 


Vs 


i> 


Tyr-1 


Gly-2 


12.6 (1.7) 


15.7 (2.2) 


4.1 (0.7) 


v 6 


LU 


Tyr-1 


Gly-2 


3.9 (0.4) 


14.4 (1.2) 


5.6 (0.5) 


V7 




Gly-2 


Gly-2 


9.1 (1.0) 


13.0 (1.4) 


4.6 (0.6) 


Vg 


i> 


Gly-2 


Gly-3 


10.6 (1.3) 


20.4 (3.1) 


7.6 (1.2) 


«9 


LU 


Gly-2 


Gly-3 


3.4 (0.2) 


16.0 (1.8) 


6.7 (0.8) 


Vio 




Gly-3 


Gly-3 


18.2 (3.2) 


31.0 (5.1) 


8.4 (1.5) 


Vll 


i> 


Gly-3 


Phe-4 


15.6 (2.9) 


52 (13) 


12 (4) 


Vl2 


LU 


Gly-3 


Phe-4 


4.4 (0.6) 


17.7 (2.7) 


7.7 (1.3) 


Vl3 


x 1 


Phe-4 


Phe-4 


3.3 (0.4) 


6.9 (1.1) 


4.4 (0.8) 


Vl4 


x 2 


Phe-4 


Phe-4 


1.7 (0.2) 


3.2 (0.4) 


3.0 (0.4) 


Vl5 


<t> 


Phe-4 


Phe-4 


8.9 (1.3) 


19.6 (3.2) 


6.3 (1.2) 


Vl6 




Phe-4 


Met-5 


4.5 (0.3) 


8.0 (0.9) 


4.4 (0.6) 


vn 


LU 


Phe-4 


Met-5 


1.8 (0.2) 


8.1 (1.2) 


5.4 (0.8) 


Vlg 


x 1 


Met-5 


Met-5 


2.7 (0.2) 


8.3 (2.5) 


6.3 (1.9) 


Vl9 


x 2 


Met-5 


Met-5 


1.9 (0.2) 


5.3 (0.4) 


4.0 (0.5) 


V20 


x s 


Met-5 


Met-5 


1.1 (0.1) 


2.7 (0.2) 


2.5 (0.2) 


V21 


x 4 


Met-5 


Met-5 


1.0 (0.1) 


1.3 (0.1) 


1.3 (0.1) 


V22 


<t> 


Met-5 


Met-5 


36 (18) 


23.8 (5.6) 


9.5 (2.4) 


V23 


<$> 


Met-5 


Met-5 


1.4 (0.2) 


1.9 (0.1) 


1.9 (0.1) 


V24 


LU 


Met-5 


Met-5 


1.0 (0.1) 


3.4 (0.2) 


3.1 (0.4) 


VO 


E 






9.0 (1.7) 


19.4 (3.1) 


6.6 (1.1) 



B. Structure 

For all our simulations we applied the method outlined 
in subsection III Bl to determine local energy minima and 
some results are summarized in Table ITVl Eq i and E'o.g 
are the lower and upper 10% quantiles of the energy and 
N C onf denotes the number of minima of the time series 
prepared for further minimization. The lowest energy 
found in this minimization process is denoted by E m i n 
and iVhits is the number of times the lowest energy con- 
figuration was hit. While the absolute values of i?o.i and 
i?o.9 vary considerably from set to set, the differences 
i?o. 9 — £-0.1 stay similar. The explanation is that the 
ASP sets differ by large additive constants to the energy. 

Again, the results of the JRF and BM solvent models 
are erratic. The BM model is entirely frozen, A^onf = 2 at 



o 




-31 -30 -29 -28 -27 -26 -25 -24 -23 
E 

FIG. 5: Local energy minima for the WE92 solvation model 
as obtained by our minimization method. 



400 K and A^onf = 1 at 300 K. Therefore, we do not give 
minimization results for BM. For JRF the N con { numbers 
are more reasonable, but still by a factor of one third and 
less smaller than the N con { numbers of each other system. 
JRF is also disregarded in the following discussion. 

Only if A^hits > 1 holds we have an indication that we 
found the GEM. Interestingly, this happens for none of 
the ASP solvent models, while it is the case for the vac- 
uum and the e(r) simulations (notably already at 400 K). 
Quite some time ago Li and Scheraga HGj developed 
a Monte Carlo minimization method and applied it to 
Met-Enkephalin in vacuum and in solvent modeled by 
OONS. While for the vacuum system their method con- 
verged consistently to the GEM, all their five runs of the 
solvent model led to different conformations with com- 
parable energies. They interpreted their results in the 
sense that Met-Enkephalin in water at 20°C is likely in 
an unfolded state for which a large ensemble of distinct 
conformations co-exist in equilibrium. A consistent sce- 
nario was later observed in NMR experiments [33|. 

Although the minimization method of Li and Scheraga 
is entirely different from ours, they essentially tested for 
valleys of attraction to the GEM at room temperature, 
quite as we do in the present paper. So, we have not 
only confirmed their old result, but find that it is also 
common to a large set of ASP models implemented in 
SMMP. Neither the method by which an ASP set was 
derived, nor whether it belongs to the fast or slow class, 
appears to matter with this respect. 

To give one example, the frequency of local energy min- 
ima of the WE92 solvation model as obtained by our 
minimization procedure from the 300 K time series is de- 
picted in Fig. N con f = 2453 minimizations are per- 
formed. Our lowest energy state is only found once and 
the same holds for the close-by low energy states. Fig. |S] 
should be compared with Fig. 2 of Ref. [jjj], where the 
frequency of the low energy minima of the vacuum sim- 
ulation is shown. There the lowest energy state relies on 
107 entries out of 1913 minimizations l34l. 
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TABLE IV: Determination of local minima: Eq.i and -Bo. 9 are the lower and upper 10% quantiles of the energies of the time 
series recorded, N con { denotes the number of configurations prepared for further minimization, E m i n (Kcal/mol) is the lowest 
energy found, and Miits is the number of times the lowest energy configuration was hit. 



Set 


T — 400 K 


T = 300 K 




Aconf 


E0.1 


E0.9 


-E'miii 


Milts 


Monf 


£0.1 


E0.9 




Miits 


VAC 


2190 


1.98 


12.26 


-12.91 


13 


1073 


-2.98 


5.73 


-12.91 


55 


e(r) 


2622 


-16.95 


-6.97 


-31.94 


8 


1312 


-21.85 


-13.17 


-31.94 


27 


OONS 


3315 


-17.83 


-9.63 


-27.69 


1 


2641 


-21.40 


-13.96 


-28.93 


1 


JRF 


448 


-317.96 


-304.98 


-328.72 


1 


365 


-323.66 


-314.24 


-332.87 


1 


WE92 


3307 


-19.88 


-11.52 


-29.44 


1 


2453 


-23.43 


-15.95 


-30.39 


1 


EM86 


2307 


8.57 


18.59 


-4.11 


1 


1191 


3.83 


12.39 


-5.47 


1 


SCH1 


2511 


5.57 


15.45 


-5.54 


1 


1147 


0.71 


9.33 


-7.52 


1 


SCH2 


3454 


-22.17 


-14.34 


-31.32 


1 


2918 


-25.31 


-18.32 


-32.71 


1 


SCH3 


2315 


8.57 


18.32 


-1.70 


1 


1229 


4.29 


12.50 


-3.29 


1 


SCH4 


2331 


8.45 


18.50 


-4.93 


1 


1108 


3.66 


12.23 


-5.16 


1 


BM 


2 


606.37 


655.16 


594.78 


1 


1 


598.37 


646.35 


590.32 


1 



VAC 300K 
VAC 400K 




-3-2-10123 
v 7 

FIG. 6: Probability density of the dihedral angle vi for the 
vacuum simulation. The arrow indicates the vacuum GEM 
value of this angle. 




-3-2-10123 
v 7 

FIG. 7: Probability density of the dihedral angle vj for the 
WE92 simulation. The arrow indicates the vacuum GEM 
value of this angle. 



In a search for structural differences of Met-Enkephalin 
in vacuum, or in the e(r) system, versus the ASP mod- 
els, we looked at the pds of the dihedral angles at 300 K. 
For all systems together there are 9 x 24 = 216 figures 
to consider. At the first look the pds of the different 
systems are amazingly similar, independently of whether 
they are from systems of the fast or slow class, from an 
ASP model, from the vacuum or from the e(r) simulation. 
A more careful investigation reveals differences, which 
appear to relate to the distinct behavior under our min- 
imization. For the dihedral angle v-j this is illustrated in 
Fig. Inland Fig.0 Its probability densities are compared 
at 300 K and 400 K. For the vacuum simulation the pds 
are depicted in Fig. HJ] and from 400 K to 300 K we ob- 
serve an increase of the peak which is located close to 
the arrow, which indicates the vacuum GEM value of V7. 



In contrast to this, the wrong peak increases in Fig. \7\ 
where the pds are shown for the WE92 solvent model. 

One may suspect that the difference between the mod- 
els of our fast and slow class is simply due to an effectively 
higher temperature for the three models of the fast class. 
To gain insight into this question, we calculate the en- 
tropies of our pds. Each pd is discretized as a histogram 
of 200 entries, py, where i = 1, ... ,24 labels the dihedral 

angles according to Table Iffll and T,f=i Pv = !• The 
entropy of the pd of a dihedral is then defined by 

200 

Si = - )^ pjj \npij (10) 

i=i 

and the total entropy of the pds of an APS model is 
S = Y^i^i- I n Fig. the thus obtained entropies are 
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400 K 














\ A 




\ fr 




300 K 















2 4 6 8 10 

itysol 

FIG. 8: Entropies of the pds of our ASP models. The models 
are labeled according Table HJ in addition itysol = for 
vacuum and itysol = 10 for the e(r) model. 



depicted for all our models. The lines between the data 
points have no other meaning but to guide the eyes. The 
dips for the JRF and the BM model show, again, that 
their configurations are essentially frozen. For the other 
we see a decrease of entropy from 400 K to 300 K, but we 
find no larger entropy for the models of the fast class than 
for the models of the slow class. Therefore, the effective 
temperature scenario is rules out. Instead, it seems that 
for the models of the fast class the solvent has some kind 
of "lubrication" effect, which accelerates the simulation. 

Strong similarities between the ASP models of the fast 
class on one side and the ASP models of the slow class 
on the other side are found for the solvation energies, the 
gyration radii and the end to end distances. 

The solvation energies E so \ yj measured during our 
simulations of ASP models are solvent-vacuum (s/v) 
transfer energies. There are structural difference between 
the typical configurations of an ASP model time series 
and the vacuum time series. Consequently, the aver- 
age s/v transfer energies are not identical with the av- 
erage vacuum- solvent (v/s) transfer energies, which are 
obtained by calculating E so \ of the solvent models on the 
configurations of the vacuum time series. The s/v as well 
as the —v/s average transfer energies are collected in Ta- 
bic The averages are taken for the canonical time 
series at 300 K and at 400 K. At 300 K averages are also 
taken for the time series minima (indicated by arrows 
in Fig. ^) and for the local minima (which are obtained 
by running the conjugate gradient minimizer on the time 
series minima). For the gyration radii R gy and the end 
to end distances i? -o the same averages are collected in 
Table (definitions and software are given in SMMP). 

For the transfer energies the over-all effect is hy- 
drophilic for the APS models of the fast class and hy- 
drophobic for the APS models of the slow class. Within 
each class the values are quite similar, despite differences 
in the interaction coefficients (see Table 1 of Ref . [l3T| ) . 
As expected the over-all transfer energies of the JRF and 



BM models are out of the reasonable range, JRF to the 
hydrophilic and BM to the hydrophobic side. Our Ta- 
ble IVII shows that we observe the extended structures 
found in previous simulations [III 

111 

and in NMR ex- 
periments [33] only for the APS models of the fast class. 



IV. SUMMARY AND CONCLUSIONS 

We have performed Met-Enkephalin simulations at 
room temperature (300 K) for the solvation models of 
Table [I] Quantitative results obtained in that way can- 
not be trusted, apparently because the methods to de- 
rive the ASPs have been quite crude. Also our simu- 
lations do not give information that would allow us to 
pick a best ASP set for the intended purpose of simu- 
lating Met-Enkephalin at 300 K. Nevertheless, we obtain 
a qualitative overview of a number of interesting con- 
sequences, which one can expect by including solvation 
effects via an ASP model in Metropolis calculations. 

Two of the ASP sets (JRF |2jJ and BM as imple- 
mented in SMMP |l3|) suffer from so large autocorrela- 
tions that for them Metropolis simulations at 300 K are 
in essence impossible. Their dihedral angles are essen- 
tially frozen. These two ASP sets are certainly erratic, 
as 300 K is a temperature at which thermodynamic fluc- 
tuations of the systems are expected (also these two sets 
perform badly at 400 K). 

The remaining nine models, seven ASP sets, e = 2 
vacuum, and an e(r) system |25j, fall into a fast and a 
slow class with respect to their integrated autocorrela- 
tion times, see Fig. 0] Vacuum simulations are in the 
slow class. This leads to the interesting feature that 
it takes less computer time to estimate physical observ- 
ables at room temperature in the fast solvation models 
OONS E3, WE92 21], and SCH2 0E2, than it takes 
for vacuum, despite the substantial increase of the com- 
puter time per sweep by a factor of about 5 for the solva- 
tion models over the vacuum system. We have no clear 
clue why some models have a fast and others a slow dy- 
namics. To derive the parameters of OONS and WE92 
vacuum- water (v/w) transfer energies were used, but for 
SCH2 it was octanol- water (o/w). Also the slow class 
features v/w as well as o/w ASP models. 

We applied the minimization procedure of Ref. 01 
in an attempt to locate the GEM for the nine systems, 
which are reasonably well-behaved under Metropolis sim- 
ulations at 300 K. The GEM is unambiguously found for 
the vacuum system and for the simulation with a distance 
dependent electrostatic permittivity. No true GEM is 
found for any of the remaining seven ASP models. This 
confirms an old result of Li and Scheraga , who con- 
cluded that at room temperature Met-Enkephalin in wa- 
ter is likely in an unfolded state. To get a better under- 
standing of this result, we studied at 300 K the dihedral 
pds in some details. At a first glance they look quite 
similar for all the models in the fast as well as in the 
slow class. Differences are found for a number of details, 
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TABLE V: Transfer energies. 



Set 


T = 300 K 


T — 400 K 




{E T ) (Local 


min.) 


(i? T )(Time series 


min.) 


{E T ) 


(E T ) 




s/v 


— v/s 


s/v 


—v/s 


s/v 


—v/s 


s/v 


—s/v 


OONS 


-24.85 


-18.71 


-24.63 


-18.53 


-24.84 


-19.07 


-24.94 


-20.19 


JRF 


-346.95 


-185.65 


-346.72 


-187.30 


-343.56 


-196.96 


-337.68 


-209.46 


WE92 


-28.84 


-17.46 


-28.52 


-17.62 


-27.87 


-18.76 


-27.47 


-20.61 


EM86 


6.19 


6.93 


6.29 


7.01 


6.45 


7.08 


6.60 


7.11 


SCH1 


2.76 


3.89 


2.86 


3.92 


2.94 


3.80 


2.97 


3.60 


SCH2 


-31.03 


-20.61 


-30.77 


-20.80 


-30.61 


-21.97 


-30.45 


-23.85 


SCH3 


4.03 


9.73 


4.28 


9.80 


4.74 


9.50 


5.15 


9.00 


SCH4 


6.08 


6.86 


6.17 


6.93 


6.32 


6.95 


6.46 


6.94 


BM 




715.29 




721.43 


581.60 


745.22 


597.61 


776.03 



TABLE VI: Gyration radii and end to end distance. 



Set 






T 


= 300 K 






T = 


400 K 




Local 


minima 


Time 


series minima 


Time 


series 


Time series 




{-Rgy} 


(-Re-e) 


(-Rgy) 


(Re-e) 


(Rgy) 


<-Re-e) 


(-Rgy) 


(Rc-c) 


VAC 


4.56 


5.67 


4.60 


5.83 


4.72 


6.83 


4.97 


8.50 


er 


4.53 


6.20 


4.57 


6.33 


4.71 


7.39 


4.99 


8.94 


OONS 


4.92 


10.30 


4.94 


10.34 


5.32 


11.71 


5.60 


12.45 


JRF 


5.75 


13.00 


5.75 


13.00 


5.78 


13.05 


5.70 


13.20 


WE92 


5.06 


12.06 


5.07 


12.09 


5.43 


13.02 


5.72 


13.47 


EM86 


4.47 


6.97 


4.48 


7.02 


4.62 


7.94 


4.86 


9.28 


SCH1 


4.51 


6.96 


4.52 


7.02 


4.68 


7.96 


4.96 


9.45 


SCH2 


5.18 


12.48 


5.20 


12.53 


5.63 


13.47 


5.86 


13.66 


SCH3 


4.54 


8.88 


4.55 


8.90 


4.72 


9.56 


4.95 


10.50 


SCH4 


4.46 


6.82 


4.47 


6.87 


4.62 


7.83 


4.87 


9.17 


BM 










4.13 


7.30 


4.21 


7.66 



which may allow to explain why the 300 K configurations 
of the ASP models behave entirely different under our 
minimization procedure than the vacuum and the e(r) 
systems. 

The central question, which remains to be settled, 
is whether ASP models will ultimately allow for accu- 
rate Metropolis simulations of biomolecules like Met- 
Enkephalin in solvent or not. In principle, this could 
be decided by determining whether ASPs exist which re- 
produce accurately mean energies of explicit solvent sim- 
ulation around a large number of fixed Met-Enkephalin 
configurations. 
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